import numpy as np
# import matplotlib.pyplot as plt
import math


def ft2(t, fz, ld, westeast):
    """
    # 参考《对应新抗震规范反应谱的功率谱模型参数研究》 周佩佩 巢 斯
    # 外包络函数
    :param t:
    :param fz:
    :param ld:
    :param westeast:
    :return:

    """
    M = 0
    R = 0
    MR1 = np.array([[(4.94, 7.740), (5.60, 7.313), (6.26, 6.896), (6.92, 6.490)],
                    [(5.60, 24.14), (6.26, 23.38), (6.92, 22.64), (7.58, 21.91)],
                    [(6.26, 52.39), (6.92, 51.02), (7.58, 49.69), (8.24, 48.39)]])
    MR2 = np.array([[(4.94, 5.337), (5.60, 5.068), (6.26, 4.824), (6.92, 4.605)],
                    [(5.60, 19.74), (6.26, 19.23), (6.92, 18.76), (7.58, 18.34)],
                    [(6.26, 44.06), (6.92, 43.08), (7.58, 42.19), (8.24, 41.38)]])
    if westeast == 0:  # 东部
        (M, R) = MR1[fz - 1, math.floor(ld - 6)]  # 7和7.5  8和8.5 落到同一组
    if westeast == 1:
        (M, R) = MR2[fz - 1, math.floor(ld - 6)]
    # 确定t1 t2 td c
    # 霍俊荣公式
    t1 = 10 ** (-1.074 + 1.005 * np.log10(R + 10))  # 上升度段终点
    ts = 10 ** (-2.268 + 0.3262 * M + 0.5815 * np.log10(R + 10))  # 平稳段长度
    # 要求：
    # 建筑：
    #   强震持续时间 至少为6s  故设置平稳长度最小为5s 基本达到要求
    # 核电：
    # 每个人工加速度时程的持续时间一般为15—30秒，其中的强信号部分（至少为最大值的25%）的持续时间应大于10秒。
    # 信号至少有6个峰值（正的或负的）超过其最大值的70%。
    t2 = max(t1 + ts, 8)
    c = 10 ** (1.941 - 0.2871 * M - 0.567 * np.log10(R + 10))

    if 0 <= t <= t1:  # 二次上升段
        kist = (t / t1) ** 2
    elif t1 < t < t2:  # 平稳段
        kist = 1
    else:
        kist = np.exp(-c * (t - t2))  # 负指数下降段

    return kist


def ft(t, cd):
    """
    # 参考《基于新《电力设施抗震设计规范》的地震动随机模型参数研究》 田 利１，高国栋１ ，盖 霞２
    # 外包络函数
    :param t: 时间
    :param cd:场地类型
    :return:
    GB 50260-2013 《电力设施抗震设计规范》:
    当采用动力时程分 进行抗震设计时，可采用实际强震记录或人工合成地震动 作为地震动输入时程。输入地震动时理不应少于 三条 至少应有一 条人士合成地震动时程。时程
     的总持续时间 不应少千 30s 其中 强震功部分不应小千 6s 计算结果宜取时程法计算结果的包络值和振型分解反应谱法计算结果的较大值
    """
    parameter_dict = {0: (0.3, 6, 0.5), 1: (0.6, 6, 0.45), 2: (0.8, 7, 0.35), 3: (1.2, 9, 0.25), 4: (1.6, 12, 0.15)}
    t1, t2, c = parameter_dict[cd]
    if 0 <= t <= t1:  # 二次上升段
        k = (t / t1) ** 2
    elif t1 < t < t2:  # 平稳段
        k = 1
    else:
        k = np.exp(-c * (t - t2))  # 负指数下降段

    return k


# # 测试
if __name__ == '__main__':
    a = np.arange(0.01, 30, 0.01)
    b1 = np.zeros_like(a)
    b2 = np.zeros_like(a)
    b3 = np.zeros_like(a)
    b4 = np.zeros_like(a)
    b5 = np.zeros_like(a)

    for i in range(len(a)):
        b1[i] = ft(a[i], 0)
        b2[i] = ft(a[i], 1)
        b3[i] = ft(a[i], 2)
        b4[i] = ft(a[i], 3)
        b5[i] = ft(a[i], 4)

    plt.subplot(511)
    plt.plot(a, b1)
    plt.subplot(512)
    plt.plot(a, b2)
    plt.subplot(513)
    plt.plot(a, b3)
    plt.subplot(514)
    plt.plot(a, b4)
    plt.subplot(515)
    plt.plot(a, b5)

    z = np.concatenate((a.reshape(len(a), 1), b3.reshape(len(b3), 1)), axis=1)
    np.savetxt('1.txt',z)



    plt.show()
